\chapter{Vibrational corrections}\label{ch:vibave}

\dalton\ provides an efficient automated procedure for calculating
rovibrationally averaged molecular $r_\alpha$ geometries, as well as an
automated procedure for calculating vibrational averages of a large
range of second-order molecular properties, for SCF and MCSCF wave
functions. In the current implementation, it is not possible to
exploit point-group symmetry, and one must ensure that the symmetry is
turned off in the calculation.

\begin{center}
\fbox{
\parbox[h][\height][l]{12cm}{
\small
\noindent
{\bf Reference literature:}
\begin{list}{}{}
\item Effective geometries: P.-O.~{\AA}strand, K.~Ruud and
P.~R.~Taylor.\newblock {\em J.Chem.Phys}, {\bf 112},\hspace{0.25em},
2655 (2000).
\item Vibrational averaged properties: K.~Ruud, P.-O.~{\AA}strand and
P.~R.~Taylor.\newblock {\em J.Chem.Phys.}, {\bf
112},\hspace{0.25em}2668, (2000). 
\item Temperature and isotope effects: K.~Ruud, J.~Lounila and
J.~Vaara. \newblock {\em J.Chem.Phys.}, to be published.
\end{list}
}}
\end{center}

\section{Effective geometries}\label{sec:effgeom}

The (ro)vibrationally averaged geometries\index{effective
geometries}\index{vibrationally averaged
geometries}\index{rovibrationally averaged geometries} can be
calculated from a knowledge of part of the cubic force field

\begin{equation}
\left<r_i\right> = r_{e,i} -
\frac{1}{4\omega_i^2}\sum_{j=1}^{3N-6}\frac{V^{\left(3\right)}_{ijj}}{\omega_j}
\end{equation}
where the summation runs over all normal modes in the molecule and
where $\omega_i$ is the harmonic frequency of normal mode $i$ and
$V^{\left(3\right)}_{ijj}$ is the cubic force field. A typical input
for determining (ro)vibrationally averaged Hartree--Fock geometries
for different water isotopomers will look like

\begin{verbatim}
**DALTON INPUT
.WALK
*WALK
.ANHARM
.DISPLACEMENT
0.001
.TEMPERATURES
 4
 0.0 300.0 500.0 1000.0
**WAVE FUNCTIONS
.HF
*SCF INPUT
.THRESH
 1.0D-10
**START
*RESPONS
.THRESH
 1.0D-5
**EACH STEP
*RESPONS
.THRESH
 1.0D-5
**PROPERTIES
.VIBANA
*RESPONS
.THRESH
 1.0D-5
*VIBANA
.ISOTOP
 3 3
 1 2 1
 1 2 2
 2 1 1
**END OF DALTON INPUT
\end{verbatim}

The calculation of (ro)vibrationally averaged geometries are invoked
be the keyword \verb|.ANHARM| in the \verb|*WALK| input module. In
this example, the full cubic force field will be determined as 
first derivatives of analytical molecular Hessians. This will be done
in Cartesian coordinates, and the calculation will therefore require
the evaluation of $6K + 1$ analytical Hessians, where $K$ is the
number of atoms in the molecules. Although expensive, it allows
 (ro)vibrational corrections to be calculated for any isotopic
species, in the above 
example for H$_2\;^{16}$O, HD$\;^{16}$O, D$_2\;^{16}$O,
H$_2\;^{18}$O. This is directed by the keyword \verb|.ISOTOP|. We
note that the most abundant isotope will always be calculated, and is
therefore not included in the list above.

We have requested that rovibrationally averaged geometries be
calculated for 5 different temperatures. By default, these geometries
will include centrifugal distortions~\cite{krjljv}. This can be
turned by using the keyword \verb|.NO CENT| in the \verb|*WALK| input
module.

By default, the numerical differentiation will use a step length of
0.0001 \bohr{}. Experience show this to be too short~\cite{poakrprtjcp112}, and we have
therefore changed this to be 0.001 \bohr{} in the example above by the
use of the keyword \verb|.DISPLACMENT| in the \verb|*WALK| input
module.

If only one (or a few) isotopic species are of interest, we can significantly
speed up the calculation of the (ro)vibrationally averaged geometries
by doing the numerical differentiation in the normal coordinates of
the isotopic species of interest. This can be requested through the
keyword \verb|.NORMAL|. The relevant part of the cubic force field is
then calculated as numerical second derivatives of analytical
gradients. We note that the suggested step length in this case should
be set to 0.0075~\cite{poakrprtjcp112}. We note that we will still need to calculate one
analytical Hessian in order to determine the normal coordinates.

The default maximum number of iterations is 20. However, \dalton\ will
automatically reset the maximum number of iterations to 6$K$+1 in case
of vibrational averaging calculations. The maximum number of
iterations can also be set explicitly by using
the keyword \verb|.MAX IT| in the \verb|**DALTON INPUT| module.

\section{Vibrational averaged
properties}\label{sec:vibavegeo}\index{vibrational
corrections}\index{zero-point vibrational
corrections}\index{temperature effects}

The change in the geometry accounts for part of the contribution to a
vibrationally averaged property, namely that due to the anharmonicity
of the potential~\cite{krpoaprtjacs123}. Although this term is important, we need to
include also the contribution from the averaging of the molecular
property over the harmonic oscillator wave function in order to get an
accurate estimate of the vibrational corrections to the molecular
property. 

At the effective geometry, this contribution to for instance the
nuclear shielding constants can be obtained from the following input
\begin{verbatim}
**DALTON INPUT
.WALK
*WALK
.VIBAVE
.DISPLACEMENT
0.05
.TEMPERATURES
 1
 300.0
**WAVE FUNCTIONS
.HF
*SCF INPUT
.THRESH
 1.0D-10
**START
.SHIELD
*LINRES
.THRESH
 1.0D-6
*RESPONS
.THRESH
 1.0D-5
**EACH STEP
.SHIELD
*LINRES
.THRESH
 1.0D-6
*RESPONS
.THRESH
 1.0D-5
**END OF DALTON INPUT
\end{verbatim}

This input will calculate the harmonic contribution to the
(ro)vibrational average to the nuclear shielding constants at 300K for
$^{17}$ODH. It is important to realize that since each isotopic
species for each temperature will have its own unique
(ro)vibrationally averaged geometry, we will have to calculate the
harmonic contribution for each temperature and each isotopic species
separately. The isotopic constitution is specified in the
\molinp\ file as described in Chapter~\ref{ch:molinp}. 

We note that we may reuse the property derivatives from a different
geometry for calculating the harmonic contribution to the vibrational
correction at the given geometry by using the keyword \Key{REUSE} in
the \Sec{WALK} module. A new force field is calculated, but the
property derivatives are assumed to remain unchanged. The
approximation has been tested and been shown to account, through the
change in the effective geometry for different temperatures, for a
very large fraction of the temperature effects on molecular
properties~\cite{krjljv}. 

This calculation will always be done in normal coordinates, and the
recommended step length is 0.05~\cite{krpoaprtjcp112}. As for the calculation of
(ro)vibrationally averaged geometries in normal coordinates, the
calculation requires the determination of one analytical Hessian in
order to determine the harmonic force field.

The default maximum number of iterations is 20. However, \dalton\ will
automatically reset the maximum number of iterations to 6$K$+1 in case
of vibrational averaging calculations. The maximum number of
iterations can also be set explicitly by using
the keyword \verb|.MAX IT| in the \verb|**DALTON INPUT| module.

It is important to understand that in some cases a property will acquire a
non-zero value only after the walk to the "effective" geometry has lowered
the symmetry of the system, an example being the dipole moment of CH3CD3:
CH3CH3 has at least $D_3$~symmetry whereas CH3CD3 at the effective
geometry has either $C_3$ or at most $C_{3v}$~symmetry, and the latter
cases permit a non-vanishing dipole moment.  A consequence of this is that
if the walk to the effective geometry does not lower the symmetry of the
Born-Oppenheimer Hamiltonian, then a property that vanishes at the
equilibrium geometry will remain zero at the effective geometry.  This is
the case for diatomic molecules.  There is only one geometry parameter,
and varying it does not change the symmetry of the electronic Hamiltonian
of the system.  A property that is nonzero at the equilibrium bond length,
such as the dipole moment of CO, or the quadrupole moment of H$_2$, can
and most likely will have a nonzero value at the effective geometry for a
particular isotopically substituted species such as HD.  But the dipole
moment of a homonuclear molecule is zero at any geometry, because of the
symmetry of the Hamiltonian and of the eigenfunctions of that Hamiltonian.
 Further, symmetry dictates that the vibrational averaging contribution
will also be zero.  Hence, for example, it is not possible to obtain the
(known) dipole moment of HD this way, because within the Born-Oppenheimer
approximation this quantity is zero.

\section{Vibrationally averaged spin--spin coupling constants}

For the calculation of vibrational corrections in indirect spin--spin
coupling constants, \dalton\ provides an alternative approach to the
calculation of vibrational corrections, in which also the full
point-group symmetry of the molecule is exploited in order to
reduce the number of displaced geometries that need to be included. An
example of such an input is :

\begin{verbatim}
**DALTON INPUT
.NMDDRV
**NMDDRV
.SYMMETRY
  C2v
.VIBANA
.DISPLA
 0.01
*PROPAV
.ANHA-P
**WAVE FUNCTIONS
.HF
*SCF INPUT
.THRESH
 1.0D-12
**EACH STEP
.SPIN-SPIN
*SPIN-S
.SELECT
  3
  1  2  3
*TRPRSP
.THRESH
1.0D-12
*LINRES
.THRESH
1.0D-12
*END OF INPUT
\end{verbatim}

We note that the input includes the full point group symmetry of the
molecule (in this case the water molecule, and thus here given as
$C_{2v}$). In this case, the vibrational corrections are evaluated at
the equilibrium geometry of the molecule, and both the harmonic and
anharmonic contributions are included. The approach is very similar to
that presented in the previous two sections, and a detailed
description of the two approaches and a comparison of the two methods
are given in Ref.~\cite{tarkr}. For information about the special input
modules used in the above input example, we refer to
Sec.~\ref{sec:nmddrv}. 
